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. Abstract. A theory of BEC interferometry in an unsymmetrical double- 

well trap has been developed for small boson numbers, based on the two-mode 



approximation. The bosons are initially in the lowest mode of a single well trap, 
which is split into a double well and then recombined. Possible fragmentations 
. into separate BEC states in each well during the splitting/recombination process 

;> are allowed for. The BEC is treated as a giant spin system, the fragmented 

. states are eigenstates of S 2 and S z . Self-consistent sets of equations for the 

amplitudes of the fragmented states and for the two single boson mode functions 
are obtained. The latter are coupled Gross-Pitaevskii equations. Interferometric 
effects may be measured via boson numbers in the first excited mode. 

o 

^ I 1 Introduction 



X 

5-S 



The realization of Bose-Einstein condensates (BEC) in cold dilute atomic gases 
■ has opened up a new area of physics research on macroscopic quantum systems, 

^ , since in a BEC at very low temperatures essentially all the bosons occupy the 

same single particle state (also referred to as modes or orbitals). Interference 
effects involving BECs were observed p], [2], and there has been considerable 
interest in various schemes for constructing high precision interferometers us- 
ing BECs [3], [4], [5]. Improvements in interferometer precision scaling as \/N 
(where N is the number of bosons) may be possible [5J. Such interferometry 
is based on the similarity between the quantum states of BECs and those for 
lasers [7J, in both cases a large number of bosons (atoms in one case, photons 
in the other) occupy a single mode, and hence BEC and laser interferometry is 
expected to be more precise than that based on single atoms or thermal light. 
The theoretical descriptions of the BEC and the laser are not quite the same 
of course. Laser light is often described in terms of coherent states (which are 
superpositions of number states), whereas in the BEC case descriptions based 
on number states are more appropriate, since superselection rules preclude su- 
perpositions of number states from being physical states [8]. In neither case 
however is the absolute phase of the laser or BEC state of any consequence for 
interferometry, indeed the idealized state of a single mode laser can be described 
by a density operator which involves a statistical mixture of coherent states with 
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all phases having equal weight, and therefore carries no more absolute phase in- 
formation than the density operator for a number state that describes a BEC. 
Absolute phase is unimportant for interferometry because interference effects 
are associated with the relative phases between two or more contributions to 
certain total amplitudes whose moduli squared determine the measured effect - 
the interferometric effects are associated with the cross terms. There are many 
forms of interferometer, but both laser and BEC interferometers just involve 
particular ways of creating such interfering amplitudes. These amplitudes may 
have different natures - in an optical Mach-Zcndcr interferometer a recombina- 
tion of two electromagnetic field amplitudes associated with splitting the EM 
field into two different spatial pathways is involved, atomic Ramsey interferom- 
eters involve combining two quantum amplitudes for a transition that can take 
place via two different quantum pathways. The interpretation of the spatial 
interference patterns seen when two independent BECs are made to overlap 
involves considering the successive detection of bosons at various spatial posi- 
tions [S], [TO], [TT], [15], [13], [H], and the interference pattern that builds up - 
which has a well-defined fringe spacing, but the absolute position of the fringes 
changes from one experiment to the next - is due to not knowing from which 
BEC any particular boson came. A well-defined relative phase is built up after 
many detections, and this is quite consistent with a fixed total boson number. 
Spatial interference effects based on successive boson detection can be described 
in terms of quantum correlation functions [15], 15], which in turn can be related 
to interfering quantum amplitudes. 

Although in principle a BEC based atom interferometer should have similar 
advantages to a laser based optical interferometer, there are effects that could 
cause problems. Firstly, unlike photons bosons interact with each other, leading 
to non-linear terms in the Hamiltonian, and this causes dephasing effects that 
could destroy the interference patterns [17], [18]. Secondly, interactions with 
the environment, single boson thermal excitations, BEC collective excitations, 
soliton or vortex formation could also cause decoherence effects. Thirdly, al- 
though it is not necessary to prepare the bosons in a coherent state to produce 
interferometric effects, nor is it necessary to develop physical elements such as 
atomic mirrors or beam splitters in exact analogy to the optical case, an actual 
process must still be designed to produce some sort of interference effect that 
is reproducible from one experiment to the next - not all interference effects 
are useful for interferometry. Fourthly, single boson detection is not as well 
developed as single photon detection, and this makes BEC interferometry more 
difficult. Fifthly, since interferometry is used for conveniently measuring other 
quantities, it is desirable that the interferometric effect should be related to the 
quantity being measured via as simple a theory as possible. 

The theory of single atom interferometers based on double well potentials is 
relatively simple [TH] , [50] , [5T] , [55] , and as interference of a BEC after splitting 
in a double well has been demonstrated [53], [53], a theory for BEC interfer- 
ometers based on such double well potentials is of some interest, and this is 
the subject of the present paper. In addition, there is a considerable theo- 
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rctical literature dealing with the behavior of BECs in double well potentials, 
describing effects such as self-trapping, Josephson oscillations, collapses and re- 
vivals of Bloch oscillations, macroscopic entanglement and so on (see [5], [25] 
for overviews). Many of these papers (see [5^1 and references therein) treat the 
BEC in a double well via various versions of a two-mode theory [57], and this 
suggests the idea of carrying out BEC interferometry in a regime where a simple 
two mode theory could be used to interpret the interferometric effects. 

The proposed BEC interferometer involves the following process. Initially 
a large number N of bosons are at very low temperature and in the same spin 
state are trapped in a single potential well in a BEC state, with all the bosons 
in the lowest mode <f>i(r). This mode is essentially symmetric. The trapping 
potential is changed from a single well into a double well and back again over 
some suitable time scale. Experimentally this might involve magnetic traps on 
an atom chip consisting of permanent magnets plus current elements, the trap 
being changed by altering a bias field. The double well potential is in general 
asymmetric and this leads to interferometric effects, such as in the probability 
at the end of the interferometric process of bosons being found in the lowest 
excited mode </>2( r )j which is essentially antisymmetric. The asymmetry in the 
trapping potential may be due to gravitational effects for example, and the idea 
behind the interferometry is to detect such asymmetry effects by measuring the 
mean number of bosons found in the excited mode. The interferometer process 
is depicted in Figure 1. 

As indicated above, the present work on double well BEC interferometry 
involves a simple theory based on the two-mode approximation. Decoherence, 
thermal, and multimode effects will be ignored and only restricted types of ex- 
citations and quantum fluctuations will be included. The theory is restricted 
to small boson numbers. Time dependent modes will be used to describe the 
adiabatic behavior, the dynamical behavior will involve amplitudes describing 
possible fragmented states of the N boson system. The system behaves like 
a giant spin system in the two-mode approximation. A variational approach 
involving spin operators will be used to determine self-consistent coupled equa- 
tions for the amplitudes and modes, the latter equations being generalizations 
of the well-known Gross-Pitaevskii equation (CPE) [28], [29] used to describe a 
single BEC. The approach is a generalization based on papers by Menotti et al 
[30] and Spekkens et al [31], both of which use variational methods. Menotti et 
al [30] however restrict the modes and state amplitudes to be Gaussian forms 
parameterized by four variational functions, and coupled self-consistent equa- 
tions are derived for these quantities. Dynamical BEC splitting, fragmentation, 
collapses and revivals are treated. Spekkens et al [31] use a variational principle 
and spin operator methods restricted to static, symmetrical potential cases to 
derive self-consistent coupled equations for state amplitudes and modes - giving 
generalized time independent Gross-Pitaevskii equations. Static BEC fragmen- 
tation is found. Cederbaum et al [35] predict fragmented excited BEC states 
in the static case using generalized time independent GPE derived using varia- 
tional methods, but restricting fragmentation to a single choice of a 50:50 split 
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between the two wells. Numerous other papers (see [55] and references therein) 
have treated BEC dynamics in a double well potential, many cither assuming 
fixed modes or that no BEC fragmentation occurs. Spin operators based on 
fixed modes have also been widely used. 

The physics of the double well BEC interferometer based on a two mode 
treatment will be discussed in section [5J The theory of the interferometer, giv- 
ing the self-consistent coupled equations for amplitudes of possible fragmented 
states and for the generalized Gross-Pitaevskii equations for the two single boson 
mode functions is presented in section [3l Considerations for numerical studies 
based on the coupled amplitude and mode equations are covered in section [H 
and the paper is summarized in section [5j Detailed quantities involved in the 
basic equations are set out in the appendix. 

2 Physics of double well BEC interferometry 

The behavior of the double well BEC interferometer involves a number of im- 
portant issues: 

1. Does the BEC fragment into two BECs (left well, right well) during the 
process? 

2. What happens to the single boson modes <j>i(r, t), 4> 2 (r, t), .as the trap po- 
tential changes? 

3. What is the essential nature of the interferometric process involved? 

4. What excited BEC states are important in the process? 

5. What effect would decoherence, quantum fluctuations, finite tempera- 
tures, .. have? 

6. How are the interferometric measurements, such as the excited boson prob- 
ability, related to asymmetry in the trapping potential? 

7. How does the interferometer sensitivity depend on the number of bosons? 

8. What is the optimum way to change the trap potential during the process? 

2.1 Fragmentation 

The possibility of the BEC fragmenting into two parts - with some bosons being 
in one mode and the rest in a second mode (see [5], [55]) - can be seen if we 
consider the energy eigenstates for N bosons in a symmetric double well poten- 
tial (see figure 2). To discuss this case we may consider two harmonic oscillator 
wells with frequency ujq separated by 2d as representing the two separate wells, 
with the actual double well having a barrier height Vb- Localized states <j> L {r) 
and <j>n{r) in each well, associated with annihilation operators Sl and or can 
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be introduced. For simplicity the extra effects due to double well asymmetry 
will be ignored at present, though of course some effects due to boson-boson 
interactions are included. 

An approximate theoretical treatment can be based on the Bose-Hubbard 
Hamiltonian - a simple model for the N boson system 

Hbh = --^{orOl + Ol^Or) 

+ + (1) 

where 

J = -2/rfr^(rr(-|^V 2 + F)^(r) (2) 
U = gJdr\Mr)\ 4 (3) 

are the tunneling and boson-boson interaction parameters. It is well-known [8] 
that there are two regimes - the Josephson regime when J 3> U and the Fock 
regime when U J. 

In the Josephson regime the ground state is given by 



E BEC = -^JN+^UN(N-1). (5) 

In this case all iV bosons are in the same delocalized state {4> L 4- <fi R )/\/2. This 
represents a ingle unfragmented condensate - the BEC phase. 

In the Fock regime the ground state is given by 



Emott = \ UN(N-2) . (7) 

In this case the two localized states cj> L and <j> R are each occupied by N/2 bosons. 
This represents a fragmented condensate - the Mott phase. 

Estimates based on harmonic oscillator wave functions 

^» = (A) 3/4 ex P (-^)exp(-(4±A (8) 

i h si/2 47rft 2 a s , n , 

ao = ( ) 9= , (9) 

mujQ m 
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gives 

J V B a d 2 

— ~ exp( 2 - ( 10 ) 

For Rb 87 with a s = 5 nm, <zq = 1 /im, cjq = 27T.58 s , Vb/TujOq — 10, we find 
J/U ~ 10~ 7 for 2d = 10 /im and J/C7 - 10+ 2 for 2d = 4 /im. Thus both 
the Fock and Josephson regimes are accessible. Hence if the interferometric 
process is adiabatic, then either a single BEC or two fragmented BECs could 
be accessed depending on the double well parameters. On the other hand if the 
process is fast, then not all adiabatic states may be accessed. For specific double 
well parameters, whether the fragmentation occurs or not will thus depend on 
the time scale of the interferometer process. The effects of asymmetry in the 
trapping potential and of more general boson-boson interactions also need to be 
taken into account, but whether fragmentation effects occur or not cannot be 
just arbitrarily assumed. 



2.2 Nature of Modes 

Since the trapping potential changes from a single well to a double well and 
back again we expect the mode functions to change during the process, and if 
the process was done very slowly the notion of time dependent mode functions 
determined via a suitable adiabatic principle is a natural one. The question 
is - what form are the time dependent mode functions likely to have? For 
simplicity the extra effects due to boson-boson interactions will be ignored at 
present, though of course effects due to double well asymmetry are included. The 
possibilities for the situation where boson-boson interactions are unimportant 
can be seen by just solving the time dependent energy eigenvalue equations [22], 
and typical results are illustrated in Figure 3. 

The situation for the single well regime is shown in Figure 3a. Here an ap- 
proximately symmetric lowest energy eigenfunction and an approximately an- 
tisymmetric lowest excited energy eigenfunction occurs, corresponding to mode 
functions at the beginning and end of the interferometer process 

In the middle of the interferometer process where a double asymmetric well 
regime occurs, two qualitatively different outcomes may occur. The two lowest 
mode functions may be approximately symmetric and antisymmetric functions 
which are delocalized over both wells. This case is shown in Figure 3b, and 
applies to situations where the asymmetry is small. On the other hand, if the 
asymmetry is larger, the two lowest mode functions are localized in different 
wells, and no longer are approximately symmetric or antisymmetric. This case 
is shown in Figure 3c. Thus, the nature of the mode functions will depend the 
trapping potential parameters, especially on the asymmetry of the double well. 
The effects of boson-boson interaction also must be taken into account, and as 
in the case of whether fragmentation effects occur or not, the form of the mode 
functions cannot be just arbitrarily assumed. 
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2.3 Interferometry Process 



Essentially, the interferometric process from t — to t = T involves an initial 
state | N, 0, 0) and a final state | N — n, n, T) representing the transfer of n 
bosons from the first mode to the second (where in general | N — m, to, t) is 
a state at time t with N — m bosons in mode (r, t) and to bosons in mode 
<p 2 {r,t)). The probability amplitude A(n,T) for the process is related to the 
transition probability via P(n,T) — \ A(n, T)| 2 and can be written in terms of 
time evolution operators U{t 2 ,ti) as 



where the transitive property of the evolution operator has been used and a 
completeness relationship involving states at time t = T/2 has been inserted. 
The last expression (|T2"j) for the transition amplitude shows it to be the sum of 
contributions at the intermediate time T/2, where m bosons have been trans- 
ferred from mode </> 1 (r, 0) to mode (j) 2 (r,T/2). Clearly, quantum interference 
in the overall transition amplitude is present, with constructive or destructive 
interference possible. In this simple exposition there are N possible quantum 
pathways present, but if the time interval between t — and t — T is divided 
into a large number of steps, the number of pathways is hugely increased. Figure 
4 illustrates the case where N = 9 and n = 1 boson is transferred into mode 
4> 2 (r,T). Here there are two quantum pathways, one where the transfer of the 
boson occurs between t = and t — T/2 and the other where it occurs between 
t = T/2 and t = T. The intermediate mode functions 4> i (r,T/2) are shown 
as localized modes, so the two intermediate states would then involve different 
numbers of bosons in the two wells. 

2.4 Excited states, decoherence, finite temperatures and 
quantum fluctuations 

Within the two-mode approximation, the basis states which can occur are lim- 
ited to fragmented states in which some of the N bosons occupy the first mode 
4> l (r, t) and the rest occupy the second mode 4> 2 (r, t) . Although superpositions of 
such states (see equations. (|3T|) . (|34)) ) can be used to describe single BEC states 
where the mode is a superposition of cj) 1 (r, t) and <j) 2 (r, t) - and such states with 
all bosons in one mode might be approximations to a collective excited state 
of the BEC - the number of collective excited states that could be described 
this way is small, yet it is known that trapped BECs have a whole spectrum of 
collective excited states (see [25], 33 J. Also, thermally excited states in which 
some of the bosons occupy further modes tfi 3 (r,t), <fi 4 (r,t), ..are also outside 
the scope of two-mode theory. Hence the two-mode theory does not allow for 
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multi-mode effects or all possible excited states that might be accessed during 
the interferometer process, especially if the initial temperature was a significant 
fraction of the BEC transition temperature. 

Decohercncc effects due to coupling with an external environment, or due to 
interactions between the BEC state and a continuum of thermally excited states, 
or due to fluctuations in the trapping potentials require treatments involving 
master equations and density operators, and this is also outside the scope of 
the pure state treatment presented here. A full theory of BEC interferometry 
taking into account excited states (collective and single particle), decoherence, 
finite temperatures, multi-mode effects and without restrictions on the boson 
number would be a worthwhile development. Such a theory could be based on 
phase space methods |34j . in which the bosonic field operator is represented by 
a stochastic space-time function, the mean value of which resembles a conden- 
sate wave function. The stochastic condensate wave function satisfies a partial 
differential equation which contains noise terms due to quantum fluctuations 
and deterministic terms resembling those in a Gross-Pitaevskii equation. Alter- 
natively, a full treatment of BEC interferometry could be based on Bogoliubov 
theory [55] , 

2.5 Interferometric measurements, sensitivity and opti- 
mum process 

Several possible interferometric effects could be measured for the double well 
BEC interferometer, including the number of bosons ending up in the excited 
mode <j> 2 (r,T) or the final spatial boson density. The objective is to find which 
responds most sensitively to the other quantities (such as gravitational fields) 
that the interferometry is intended to measure, and this can only be determined 
via numerical studies of the operation of the interferometer. Such studies will 
include varying the parameters describing the process, such as the time scales, 
barrier heights, separation of the double wells, boson numbers and so on, to 
maximize the interferometric effects. 



3 Theory 

In terms of bosonic field operators ^(r), 1 !^ (r) the Hamiltonian is given by 

H = fdr (— V$ f • W + $ t y§ + ^§t$t§§\ ( 13 ) 
J \2m 2 J 

The first term represents the kinetic energy of the bosons each of which has mass 
m, the second term involves the time-dependent trapping potential V(r,t) and 
the third term allows for the two-body interaction between the bosons in the 
usual zero-range approximation. The coupling constant g is determined from 
the scattering length o s via g = \-na s T? jm. Since a single component BEC is 
involved only one pair of field operators is required. 
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The field operators satisfy the usual bosonic commutation rules 
^(r),*^!-')] = S(r-r') 



(14) 



Time dependent single boson mode functions <^j(r,t) will be used, chosen to 
be orthogonal and normalized at all times. 



/dr^*(r,t)(^(r,f) =S t 



(15) 



The conditions in equation (|15|) for each time t will act as constraints in the 
variational method used to obtain equations for the two mode functions. 

The field operators are expanded in terms of the mode functions, which in- 
troduces the mode annihilation Ci(t) and creation operators Ci (t) as the time 
dependent operator expansion coefficients, the mode functions carrying all the 
position dependence. The creation and annihilation operators satisfy the stan- 
dard bosonic commutation rules at all times. 



*(r)= £ 5(*)<fc(r,t) $ f (r)= £ c^)0*(r,i) 



i=l,2 



i=l,2 



= 1,2,..) 



(16) 



(17) 



In the two-mode approximation only two terms are included in the expansions 
for the field operators. 

The boson number operator N is defined by a space integral involving the 
field operators and may be also expressed as a sum involving mode annihilation 
and creation operators. Thus: 



N = /dr* (r)*(r) 



(18) 
(19) 



The boson number is a conserved quantity and only state vectors with a single 
boson number N will be considered here. For convenience N will be even. 
In a two-mode theory it is convenient to introduce spin operators defined by 

S x = (&ci + ci'c 2 )/2 

S y = (cVci - ci f C2)/2i (20) 
S z = {c 2 ] c 2 - ci f ci)/2 

The spin operators S a satisfy the standard commutation rules for angular mo- 
mentum operators 



ie a0y S-y (a,/3,j = x,y,z), 



(21) 



9 



and the square of the angular momentum (S) can be related to the boson 
number operator. Thus: 

(Sf = US a ) 2 (22) 

- |(f + 1 ) 

Clearly the angular momentum squared is a conserved quantity. 
A set of states for the N boson system can be defined by 

\k) = ^ r {k = -N 2,-N 2+1,. .,+N 2) (24) 

Kf -fc)!] 1 [(f + k)l]$ 

In general this represents a state with (^ — k) bosons in mode 4>i(r, i) and 
(5r + k) bosons in mode 4> 2 {r, t). Such a state is a fragmented state of the 
N boson system, involving two BECs not just one. These states will be used 
as orthogonal, normalized basis states for representing a general state of the 
bosonic system during the interferometer process. For the cases where k = 
±iV/2 the N bosons are all in the same mode, so that an unfragmented single 
BEC is represented. Thus with k = —N/2 we have 



2 / [JV!]j 



0) . (25) 



This state is a single unfragmented BEC with all bosons in mode (f> 1 (r, t) . 

The N boson system behaves like a giant spin system in the two-mode ap- 
proximation. The basis states | k) are simultaneous eigenstates of (S) 2 and S z 
with eigenvalues y(t + ^ an< ^ ^' Thus: 

(S) 2 \k) = ^ + l)\k) (26) 
S z \k) = k\k). (27) 

Hence j = ^ is the spin angular momentum quantum number, and k is the 
spin magnetic quantum number, with (— y < k < y). Thus the boson number 
N and the quantity k that specifies the fragmentation of the BEC between 
the two modes have a physical interpretation in terms of angular momentum 
theory. Since boson numbers may be ~ 10 s the spin system is on a macroscopic 
scale. To emphasize the spin character of the basis states we can introduce the 
notation 

' N 



k) 



(28) 



The methods of angular momentum theory can be utilized by first writing 
the Hamiltonian in terms of spin operators using equations (fTo| . (|2U|) . and 
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its matrix elements calculated using angular momentum theory from previous 
expressions plus 



S± 



N 

r 1 



r N ,N 



fc(fc±l)}5 



N 



,k±l 



S± — S x i iSy. 



(29) 
(30) 



The quantum state | of the N boson system during the interferometer 

process will be written as a superposition of the fragmented states | fc) , where 
the amplitude for this fragmented state is bk(t). 

N 

\Ht))= k i N h(t)\k). (31) 

Normalization of the state vector requires that the amplitudes satisfy the con- 
dition 

N 

t Mt)? = l, (32) 
fe=-# 



which represents conservation of probability. The condition in equation (|32[) 
for each time t will act as constraints in the variational method used to obtain 
equations for the amplitudes. The initial condition involves having a single BEC 
with all bosons in mode ^> 1 (r, 0), thus: 



$(0)) 



The form of the state vector given in equation (|3"Tj) involves a physical as- 
sumption in that only the two mode fragmented states are included in the 
quantum superposition. This amounts to ignoring other possible states for the 
bosonic system, such as where bosons occupy more than two modes or where 
collective excited states such as breathing modes are involved. Further develop- 
ment of the theory to allow for the presence such other states may be required 
if the present simple approach proves inadequate. 

It should be noted that as well as allowing for the possibility of fragmentation 
of the BEC into two modes, the state vector in equation (|3"Tj) is also consistent 
with the situation where all N bosons are in a single mode of the form 

4> 1 = cos9 exp(-i^x) 4>i + sin 6< cxp(+i^x) <i> 2 , 

where 9 determines the relative contributions from the original modes (f> l and 
(f> 2 , and where x is a phase variable. In this case the amplitudes bk are related 
to binomial coefficients and are given by 

(cos6»)- _fc (sin6»)- +fe exp(-ifcx). (34) 



bk 



(f-*0'(f 



k)\ 



11 



This situation amounts to replacing the two mode functions 



by 



(where (f> 2 = —sm8 exp(—i^x)4>i + cos# exp(+i^x) <p 2 ). The state vector is 
then given by an expression analogous to equation (jiM)) with k = —N/2, but 
with the original creation operators c\\ c 2 replaced by new creation operators 
associated with the new modes cf> 1 , <fi 2 . If it turns out that the BEC does not 
fragment then the solutions for the amplitudes b k will be in a form given by 
equation (|34[) . Such states with all bosons in one mode might approximately 
represent a collective excited state of the BEC. 

The amplitudes b k (t) and the mode functions (j>i(r,t) can then be related to 
the various types of interferometer measurement. For example, the number of 
bosons in the mode (f> 2 (r, t) is given by 



N 2 



<S>(t)\cl(t)c 2 (t)Mt) 



N 9 

= T +£*N 

1 k 



(35) 
(36) 



The time dependence is left understood in the result. Measurement of N 2 at end 
of the process depends on the asymmetry and exhibits interferometric effects 
because the probability amplitude at the end of the process for fragmented 
states with k ^ —N/2 in which there are bosons in the mode 4> 2 (r,t) will 
contain contributions from many quantum pathways. Interferometric effects of 
the spatial type can be described in terms of quantum correlation functions [15] , 
[16] . For example, the first order correlation function is given by 



G (1) (r,r',<) = (<$>(t)\¥(r)^(r')\$(t)) (37) 
= E b k *b k {^(r)*<Mr') (y - fc ) + Mr)* Mr') (y + k 

+ E b k *b k+1 J <Mr)> 2 (r')W (f - fe) (j + k + lj I 



+ E b k *b k _, J <t> 2 {vTUv')J (^ +k ^(^-k + l 



(38) 



where in the result the time dependence is left understood. More complex 
expressions are involved for the second order correlation function. The presence 
of spatial interferometric patterns and the existence of long range order in BECs 
can be determined from such correlation functions. 

The equations governing the amplitudes b k (t) are obtained from a variational 
principle based on the dynamical action Sd yn - This quantity is a functional of 
quantum state \&(t)) and is defined by 



S dyn =Jdt ({(0 t $| $) - ($| dt$)}/ 2i - ($| H |$) / h 



(39) 



12 



The Principle of Least Action involves the minimization of the action Sd yn 
for arbitrary variations of the state vector and this results in |$(t)) satisfying 
the time-dependent Schrodinger equation (TDSE). The variations of the state 
vector are subject to the constraint that it remains normalized to unity. This 
variational principle may be regarded as the fundamental principle of quantum 
dynamics, so its application to a specific case such as the BEC interferometry 
process is on firm ground. In the present situation the state vector is restricted in 
its possible variations to remaining in the form given in equation (|3ip (though 
remaining normalized to unity), and hence does not itself satisfy the TDSE. 
What is obtained is a state vector which is an approximate solution to the 
TDSE, and it turns out that the amplitudes bk(t) involved in the form for 
the state vector could also be obtained by just assuming that |$(f)) satisfied 
the TDSE. The present variational approach has been applied in many other 
quantum physics problems - the derivation of the time-dependent Hartree-Fock 
equations for electrons in an atom being one example. It has already been 
applied to BEC problems by Menotti et al [301, who described the amplitudes 
via a Gaussian function with two variational parameters. 

For fixed modes <^(V, t) the action Sdyn is a functional of the amplitudes 
bk(t). The normalization constraint in equation (|32p for time r may be written 
in terms of the functional F T [bk, b* k ] 1 which is required to equal unity. Thus 

F T [b k ,bt\ = fdt £ &?(*)&,(*)*(* -t) = 1. (40) 
i 

The action Sd yn is minimized for arbitrary variation of the amplitudes subject to 
the normalization constraints, which are taken into account with Lagrange mul- 
tipliers X(t)/H. In applying the Principle of Least Action, the functional deriva- 
tives of the action Sd yn plus the integral of the constraints F T each weighted 
with Lagrange multipliers X(r)/h are equated to zero. Thus we have: 

-^AS dyn [b k ,bt) = A.AS dyn [b k ,bl}=0 (41) 

AS dyn [b k) b* k ] = S dyn [b k ,b%] + JdT^F T [b k ,bl] (42) 

It turns out that the Lagrange multiplier A(r) associated with the normalization 
constraint can be transformed away and need not appear in the equations for 
the amplitudes. The key equations for the amplitudes b k (t) are given below in 
equation (|47|) . 

The equations governing the mode functions (j>i(r, t) are also obtained from a 
variational principle, but now based on the adiabatic action S a dia- This quantity 
is a functional of quantum state |$(t)) which is defined by 

S adia =Jdt(-($\H\$)/h) (43) 

This second Principle of Least Action involves the minimization of the action 
Sadia for arbitrary variations of the state vector, and this results in |$(t)) satis- 
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fying the time-independent Schrodinger (or energy eigenvalue) equation (TISE). 
The variations of the state vector are subject to the constraint that it remains 
normalized to unity. This variational principle may be regarded as the fun- 
damental principle for determining energy eigenstates, so its application to a 
specific case such as the BEC interferometry process is on firm ground. As 
before, the state vector is restricted in its possible variations (though remain- 
ing normalized to unity) to remaining in the form given by equation (|3ip . and 
hence does not itself satisfy the TISE. What is obtained is a state vector which 
is an approximate solution to the TISE. However, the time-dependent mode 
functions that are obtained from the variational principle can not be obtained 
just by substituting for |$(£)) in an energy eigenvalue equation. This varia- 
tional approach has been applied in many other quantum physics problems - 
the derivation of the standard time-independent Gross-Pitaevskii equation for a 
single BEC being one example. It has already been applied to other BEC prob- 
lems involving symmetrical double well potentials by Spekkens et al |31j . The 
application of the Least Action Principle to the adiabatic action to determine 
the mode functions and to the dynamic action to determine the amplitudes is 
designed to produce mode functions that would apply if the trapping potential 
were to change adiabatically, and to generate amplitudes that describe dynami- 
cal behavior in which the bosonic system may involve changing superpositions of 
different fragmented states. However, as will be seen below, the mode functions 
also reflect the possible way the BEC could fragment, with the more important 
fragmentation possibilities having greater influence in determining the mode 
functions. This is more realistic than determining mode functions based on 
some a priori assumption about fragmentation. 

For fixed amplitudes bk(t) the action S a dia is a functional of modes <f>i(r, t). 
The orthogonality and normalization constraints in equation (|15p for time r 
may be written in terms of the functionals Gf[(j) i ,(f>*], which are required to 
equal 8u- Thus 

G?[<j> it <f>*} =fdtj dr#(r,i) Mr, t) S(t - r) = 6 kl (44) 

The action S a dia is minimized for arbitrary variation of the modes subject to the 
orthonormality constraints. The functional derivatives of the action S a dia plus 
the sum, integral of the constraints G™' each weighted with Lagrange multipliers 
Nfi k i(T)/H are equated to zero. Thus we have: 



6 S 

j— AS a dia[M<f>*} = T—ASadia[<j>i,<f>*}=0 (45) 
AS a dia[<Pi,(p*] = Sadia[<j>i, <Pi] + 

+ J2jdr^tlGf[ <t > i , ( f>;] (46) 

kl 

The Lagrange multipliers associated with the mode orthonormalization con- 
straints form a Hermitian matrix of generalized chemical potentials fi^ (t) . The 
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key equations obtained for the modes (\> i (r, t) are coupled generalized Gross- 
Pitaevskii equations and are given below as equation (|48|) . These equations 
are time-independent in that no time differentiation of the mode functions is 
involved, but they are time-dependent because the mode functions are time- 
dependent due to the presence of the time-dependent trapping potential V(r,t). 

The coupled amplitude equations obtained are 
rih 

= E( H ki - nU M )bi (fc = -N/2, .., N/2). (47) 

These N+l equations (|4"T|) describe the system dynamics as it evolves amongst 
the possible fragmented states. The equations are similar to the standard ampli- 
tude equations obtained from matrix mechanics. In these equations the matrix 
elements H^i, Uki depend on the mode functions <^(r, t). Detailed expressions 
for Hki , Uki are given in Appendix [6] The matrix elements Hki are in fact the 
matrix elements of the Hamiltonian H in equation (|13p between the fragmented 
states | k) , 1 1) . The matrix elements Uki are elements of the so-called rotation 
matrix, and allow for the time dependence of the mode functions. 

The coupled equations obtained for the two modes are 

j j ^=x,y,z 

+SE^mn^J B (*=1,2). (48) 

jmn 

These two equations (|48|) describe the adiabatic behavior of the two modes. 
The equations are coupled generalized Gross-Pitaevskii equations, rather than 
the usual single mode Gross-Pitaevskii equation [28] , [29] . The coefficients Xy , 
Yij mn depend quadratically on the amplitudes bk(t). The Xy are ~ N, and the 
Yijmn are ~ N 2 . Detailed expressions for Xy, Yy mn are given in Appendix 
[S] The quantities fi^ form a 2x2 Hermitian matrix to be referred to as the 
chemical potential matrix. Together the combined set of equations for the am- 
plitudes and modes form a self-consistent set - neither the amplitude equations 
nor the generalized Gross-Pitaevskii equations can be solved independently of 
the other. This self-consistent feature is absent from most other treatments of 
BEC dynamics - the fragmentation behavior is often studied assuming that the 
modes are known in advance and considered fixed, whilst the mode functions 
are often calculated assuming some specific fragmentation, such as having half 
the bosons in each well. In the present work, the generalized Gross-Pitaevskii 
equations reflect the relative importance of all the possible fragmentations of 
the N bosons into the two modes. 

The energy E of the bosonic system can also be expressed in terms of the 
mode functions ^(r, t) and amplitudes We find that 



15 



<$(t)|ff|$(t)> (49) 



ij £TYl ^—x,y.z 

+ ~ E Yij mn / * 4>*i 4>*j <P m <f>n- ( 50 ) 
^ ijmn 

As can be seen, the energy also depends on coefficients Xy, Yij mn . 

The chemical potential /i is defined as the derivative of the energy with 
respect to the boson number, and roughly gives the change in energy if one 
boson is added to the system. By writing = xh 'N + 0(N°) and Yij mn 

= y^mnN 2 + 0(N 1 ) an expression for the chemical potential can be obtained 
using equations (f50|) . (|48|) . Thus we have 



dE 
~dN 

E^ + O(iV ). (52) 



(51) 



This result shows that the form a generalized chemical potential matrix, the 
trace of which is the chemical potential. 

The initial conditions for the amplitudes in the case where all the bosons are 
in mode <p l will be 

h(0) = 5 k> _N. (53) 

In this case only non-zero coefficients are 

X n (0) = N Y nu {0) = N(N-l), (54) 

and all the chemical potential matrix elements all zero except for /i n . We find 
that the mode function <^) 1 (r, 0) at time zero will then satisfy a single Gross- 
Pitaevskii equation of the form 

Mu ^ = -A- E Btfa + Vfa + giN-l) |<K| 2 K (55) 

This result is the expected one for the case where all bosons are in mode 4>i- 
The other mode function 4> 2 (r,0) is chosen by orthogonality. 

The regime of validity for the present two-mode theory is determined using 
the criteria that the mean field energy N g \4>\ 2 is small compared to trap phonon 
energy fnoo [36] . and the temperature T is much smaller than the transition 
temperature T c . Applying these criteria lead to conditions on the boson number 
N and the temperature T 

N « ^ (56) 
a s 

T « 0.94^/3^, (5?) 
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where clq = \f (K/2mu>o) is the harmonic oscillator vibrational amplitude. For 
Rb 87 with a s = 5 nm, a = 1 fim, uj = 2vr.58 s"\ find N < 2.10 2 and T < 
15.4 nK. Evidently the boson system can not be too large, nevertheless these 
conditions are realizable. Boson detection would be facilitated using metastable 
He 4 to form the BEC. 

4 Numerical Studies 

Numerical solutions for the amplitude and generalized Gross-Pitaevski equa- 
tions (|47|) . ()48|) involve representing the amplitudes on a time grid and the 
mode functions on a space-time grid. The calculations would be facilitated by 
introducing dimensionless units for space and time based on harmonic oscillator 
units. 

If there are Nt time points and Nsx ,Nsy ,Nsz space points for each of the 
three space dimensions respectively, then the amplitudes and the mode functions 
will require (N + 1)Nt and 2NtNsx-Nsy -Nsz complex values respectively - in 
all Nt(N + 1 + 2Nsx-Nsy-Nsz) values. The chemical potential matrix would 
also require another ANt values. Initial studies will be for the case where the 
splitting is essentially in one direction (Z), with the system tightly trapped in 
the two transverse (X, Y) directions. In this case it may be sufficient to take 
N sx = N SY = 10 and N S z = 10 3 . With N T =!0 3 systems with up to about 
N = 10 5 bosons would require about 3x10 s values if all time or space-time 
values for amplitudes, mode functions, chemical potentials were to be stored in 
the computer. 

Two possible approaches to carrying out the numerical studies are as fol- 
lows. Both involve an iterative process. These may be referred to as: (a) Time 
evolution method (b) Matrix method 

4.0.1 Time evolution method of solution 

First Step: 

1. Assume the amplitudes bk(t), the mode functions </>, ; (r, t) and an initial 
choice of their time derivatives 9j^(r, t) are known at time t 

2. Calculate the spatial derivatives of the mode functions via 

S^(r, t) ~ (^(r + Ar M , f) - 0.(r, t))/Ar„ (58) 

3. Calculate the Hki(t) from (jTTj) using equations ([63]) . ([64]) for Wy(r, t)and 
Vijmn(r,t) and calculate Uki{t) from (f68|) using (f65|) for Ty(r, t) 

4. Use the approximation for small At 

b k {t + At) ~ b k (t) + E(H k i(t) - hU M (t))bi(t) (59) 
in i 
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together with applying the normalization requirement (|32[) to determine 
the amplitudes bk (t + At) at time t + At 

Second Step: 

1. Calculate the Xij(t + At) and Yij mn (t + At) at time t + At from equations 

2. Solve the generalized GPE (|48|) for the mode functions <j>i(r,t + At) at 
time f + At 

Third Step: 

1. Improve the values of the time derivatives 9t<^j(r, t) at time t via the 
expression 

d t <t>i(r,t)~ (^(r,t + Ai)-^(r,t))/At (60) 

2. With the new c? t </>j(r, i) at time t go back to the first step and iterate the 
process until these time derivatives converge 

3. The final dt<f>i(v, t) may then be used as the initial choice for 9t^(r, t+At) 
at time t + At 

Fourth Step: 

1. As the amplitudes bk(t + At), the mode functions 4>i(r,t + At) and an 
initial choice of their time derivatives dt4>i{Y,t + At) are now known at 
time t + At we can go back to the first step and repeat the process to 
obtain the results at time t + 2At 

2. The process continues for further time points t + 3 At, t + 4 At, t + 5 At, .. 
Fifth Step: 

1. The process begins with t = using the initial amplitudes 6^(0) given by 
(|53|) and mode functions ^(r, 0) obtained from (|55|) and orthogonality. 
The initial choice of time derivatives at t = may be assumed to be zero, 
as the process will correct this initial arbitrary choice. 

The advantage of the time evolution method is that the values for the 
amplitudes, mode functions, their spatial and time derivatives and the chem- 
ical potentials need only be retained at two times t and t + At, thus only 
2(7V + 5 + 10Nsx-Nsy -Nsz) simultaneous values would be stored. If we take 
Nsx = N S y = 10 and N S z = 10 3 , then systems with up to about N = 10 5 
bosons would require about 2xl0 6 values to be simultaneously stored in the 
computer. 
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4.0.2 Matrix method of solution 

First Step: 

1. Assume a solution for the amplitudes bk as functions of time 

2. Calculate the Xij and Yij mn as functions of time 

3. Solve the generalized GPE (|48[) for the mode functions 4> i as space-time 
functions via non-linear matrix methods 

Second Step: 

1. Using equations (|55|) . (|60[) to obtain the spatial and time derivatives, cal- 
culate the Hki and Uki as functions of time 

2. Solve the amplitude equations (|4"7|) for the amplitudes bk as functions of 
time via matrix methods. 

Third Step: 

1. Repeat the process until the solutions for the mode functions and ampli- 
tudes converge. 

This approach represents the space-time values and time values of the mode 
functions and amplitudes in a column vector and then the non-linear equations 
for this vector obtained from equations (|47|) . (|48|) are solved via matrix meth- 
ods. Here the values for the amplitudes, mode functions, their spatial and time 
derivatives and the chemical potentials need only be retained at all times, which 
as we have seen would require about 3x10 s values for systems with up to about 
N = 10 5 bosons. 

5 Summary 

Using the two-mode approximation and treating the N bosons as a giant spin 
system, a theory of BEC interferometry has been developed by applying the 
Principle of Least Action to a variational form for the quantum state which 
allows for the possibility that the BEC fragments into two, as well as for the 
outcome where only a single BEC ever occurs. The amplitudes for the possible 
fragmented states describe the dynamics and are determined from the dynamic 
action. The two spatial mode functions describe the adiabatic behavior and are 
obtained from the adiabatic action. 

Self-consistent coupled equations have been obtained for the state amplitudes 
and the modes, the former being in the form of standard matrix mechanics 
equations, the latter equations being a generalization of the time independent 
Gross-Pitaevskii equations and which involve generalized chemical potentials. 
The self-consistent feature is that the mode functions are needed to determine 
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the Hamiltonian and rotation matrices that appear in the amplitude equations, 
whilst the amplitudes for possible fragmented states determine coefficients that 
appear in the generalized Gross-Pitaevskii equations for the modes. Unlike 
previous work, the mode equations reflect the relative importance of all the 
possible divisions or fragmentations of the bosons into two modes. 

Numerical studies of these equations are planned, aimed at applications in 
future BEC interferometry experiments at Swinburne University of Technology 
involving a double well interferometer based on atom chips. Two approaches 
for carrying out these numerical studies have been outlined. 

6 Appendix - Expressions for quantities in am- 
plitude and mode equations 

In the two-mode approximation the N boson system behaves like a giant spin 
system with spin quantum number j = N/2 and which can be described via 
angular momentum eigenstates | fc), where k — —N/2, .., +N/2 is a magnetic 
quantum number which describes fragmented states of the bosonic system with 
— k) bosons in mode (f>i(r, t) and + k) bosons in mode <j> 2 (r,t). It is 
therefore not surprising that the basic equations will involve expressions arising 
from angular momentum theory. These are the quantities and YJy mn which 
are defined as 
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= {(y-0(-j + ^«=(y+*)*« (61) 



\rll 11 




AT 

2-*-l)*« 


T/ Z»i 

x kl 






v 12 12 
x kl 
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(62) 



The Hamiltonian and rotation matrix elements Hki and t/fc; that occur in the 
amplitude equations [)47p involve spatial integrals involving the mode functions 
<fi l and They are therefore functionals of the mode functions. The expres- 
sions depend also on the^ spatial and time derivatives of the mode functions 
through the quantities Wy(r,i), Vy TOn (r,i) and Ty(r,i), where (i,j,m,n — 
1,2), and which are defined by 

%(r,t) = £ d^d^+tfVh (63) 



9 ?77 Z " 

9 i* x* 

2 C 
1 

2i' 



Kjmn(r,£) = 7;<f>i<l>j<P m <Pn ( 64 ) 
T y (r,t) = ^(dttfh-fidth) (65) 



The rotation matrix elements £7*,; (— y < k, I < + y) are given by 

Uh = ±-.[(d t (k\)\l)-(k\(d t \l))} = U? k (66) 
= f drUufatfAhAti). (67) 
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In the expression (|67[) for the rotation matrix the quantity Uu is 

Uki = J2 X kiTir (68) 

ij 

The result involves the angular momentum theory quantities XS . Thus for the 
rotation matrix, space integrals of the mode functions and their time derivatives 
are involved. 

The Hamiltonian matrix elements Hki (— -y < k, I < +-^) are given by 

H k i = (k\H\l)=H? k (69) 
= fdrHufcJld^dnti). (70) 

In the expression (|70p for the Hamiltonian matrix the quantity fffc; is a Hamil- 
tonian density and is given by 

H kl =J2 X klWv+E Y M " m % mn- (71) 

ij ijmn 

This result involves the angular momentum theory quantities XA and Y^j" 171 . 
Thus for the Hamiltonian matrix, space integrals of the mode functions and 
their spatial derivatives are involved. 

The coefficients Xy and Yij mn (i,j,m,n — 1, 2) that occur in the generalized 
Gross-Pitaevskii equations ([48]) for the mode functions are quadratic functions 
of the amplitudes bk (— y <k,l < +^) 

^ = E^4i^=^~^ (72) 

n,™ = E6i^ wm 6i=^«~JV 3 (73) 
lb,i 

Note the Hermitian properties of these quantities and the N dependence of their 
order of magnitude. 

7 Figure captions 

Figure 1. The interferometer process. A trapping potential (shown in red) is 
changed from a single well into an asymmetric double well and back to a single 
well again. Initially all the bosons (shown as squares) are in the symmetric 
lowest mode of the single well, at the end of the process some bosons are in the 
antisymmetric first excited mode of the single well. Mode functions are depicted 
in pink and blue, and possible changes to the mode functions during the double 
well intermediate stage are shown. 

Figure 2. Bosons in a symmetric double well trap showing possible frag- 
mentation effects. For low barrier heights and small inter-well separation (as in 
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(a)) a single unfragmcntcd BEC occurs, with all bosons in the symmetric mode 
delocalized between the two wells (Josephson phase). For the opposite situation 
(as in (b)) the BEC fragments into two, with half the bosons in localized modes 
in each well (Mott phase). Trap asymmetry is ignored. 

Figure 3. Mode functions in asymmetric trapping potentials showing lo- 
calization and derealization effects in the double well regime. For the single 
well regime (a) the symmetric and antisymmetric two lowest modes are shown. 
For the double well regime with small asymmetry (b) two delocalized modes are 
shown, one approximately symmetric the other approximately antisymmetric. 
For the double well regime with large asymmetry (c) two localized modes are 
shown, each localized in a different well. Boson-boson interactions are ignored. 

Figure 4. BEC interferometry as a quantum interference process. The case 
with N — 9 bosons initially in mode (f) 1 (r, 0) and n = 1 bosons finally transferred 
to mode <f> 2 (r,T) is shown. Two quantum pathways are present depending on 
whether the transfer occurs between t = and t = T/2 or between t = T/2 and 
t = T. 
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